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ABSTRACT 

The formation and evolution of a wide class of astrophysical objects is governed by turbulent, magnetized 
accretion disks. Understanding their secular dynamics is of primary importance. Apart from enabling mass ac¬ 
cretion via the transport of angular momentum, the turbulence affects the long-term evolution of the embedded 
magnetic flux, which in turn regulates the efficiency of the transport. In this paper, we take a comprehensive 
next step towards an effective mean-field model for turbulent astrophysical disks by systematically studying the 
key properties of magnetorotational turbulence in vertically-stratified, isothermal shearing boxes. This allows 
us to infer emergent properties of the ensuing chaotic flow as a function of the shear parameter as well as the 
amount of net-vertical flux. Using the test-field method, we furthermore characterize the mean-field dynamo 
coefficients that describe the long-term evolution of large-scale fields. We simultaneously infer the vertical 
shape and the spectral scale dependence of these closure coefficients, with the latter describing non-local con¬ 
tributions to the turbulent electromotive force. Based on this, we infer a scale-separation ratio of about ten 
for the large-scale dynamo. We finally determine scaling properties of the mean-field dynamo coefficients. 
The relevant component of the dynamo a effect is found to scale linearly with the shear rate, as is the corre¬ 
sponding turbulent diffusion, //. Together, these scalings allow us to predict, in a quantitative manner, the cycle 
period of the well-known butterfly diagram. This lends new support to the importance of the ail mechanism in 
determining the evolution of large-scale magnetic fields in turbulent accretion disks. 

Subject headings: accretion disks - magnetohydrodynamics (MHD) - methods: numerical 


1. INTRODUCTION 

Accretion disks are of central importance to the formation 
and evolution of a wide range of astrophysical objects. If 
sufficiently ionized, such disks are subject to the magnetoro¬ 
tational instability (MRI), driving turbu lence within the disk 
body. Since the pioneering work of|Balbus & Hawley|<|1991|>, 
numerical simulations have played a key role in establishing 
our current kno wledge about the non-lin early saturated, tur¬ 
bulent state (see |Balbus & Haw ley {2003 for a review) - yet 
there remain substantial open questions regarding the satura¬ 
tion amplitude of the turbulent stresses, and their scaling with 
various input parameters. 

Ideally, one would strive to simulate the accretion disk in 
its entirety, but the numerical requirements to treat turbulent 
magnetized disks globally are immense even by today’s stan¬ 
dards (see Qgilvie|2003[ for an estimate). Both the grid reso¬ 
lution ( [Hawley et al.|2013| ) and the domain size (Flock et al. 
2012)1 have to be sufficient for obtaining converged results, 
and existing global models (Hawley] |2000) | Fromang & Nel- 

Beckwith 


|son|2009HDzyurkevich et al.|2010) Beckwith et al.|2011[ ) all 

have to make compromises of one sort or another. 

The need for global simulations (necessary to properly de¬ 
scribe the disk wind and jet) on one hand, and the diffi¬ 
culty to resolve all relevant length scales (necessary to ob¬ 
tain converged results for the turbulence) on the other hand, 
strongly promotes the development of a so-called sub-grid- 
scale model (in the vein of the large-eddy simulations pio¬ 
neered by Smagorinsky |1963) . Such an approach can be re¬ 
garded as mandatory when integration over long dynamical 
timescales is required, for instance to follow the secular evo¬ 
lution of the accretion disk, or the stab ility and symmetry of a 
protostellar jet (Fendt & Sheikhnez ami||2013] >. 
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Here we advocate a strategy where adequately resolved 
local shearing-box models with varying shear rate, q, and 
net-vertical flux (NVF) are employed to derive an effec¬ 
tive parametrization of MRI-generated turbulence with the 
prospect of later inserting this into global simulations as a 
sub-grid-scale model. Early phenomenological descriptions 
for the turbulent stresses (jKato & Yoshizawa][l995} JOgilvie 
2003) included the effects of the background shear but ig- 
ed a 


ig¬ 
nored any large-scale fields. These phenomenological ap¬ 
proaches typically involve a certain number of unc onstrained 
order-unity c oefficients. Liljestrom et al. (2009) and Snellman 


et al.| (|2012|) have made first attempts to restrict various clo- 
sure coefficients in heuristic approaches, such as the Ogilvie| 


(2003) stress model, from direct simulations, pioneering an 
important future route. Penna et al. ( 201 3| ) have tested such 
shear-stress relations against global relativistic simulations, 
and derive a prescription where the turbulent stress depends 
on radius. Starting from the early phenomenological models, 
|Pessah, Chan, & Psaltis| ( 2006a] ) included for the first time a 
driving term akin to the linear MRI assuming a constant net- 
vertical flux. Notably, this simple model already produces a 
scaling of the Maxwell-to-Reynolds stress ratio according to 
(4 - q) /q as a function of the dimensionless shear parameter 
( [Pessah, Chan, & Psaltis|2006b|), as well as the tot al stress as 
a function of shear ( [Pessah, Chan, & Psaltis||2008| >. The for¬ 
mer relation had previo usly been found in unstratified MRI 
simulations by Hawley, Balbus, & Winters j\999) . This cor¬ 
respondence may indicate that the linear modes remain impor- 
tant even when affected by turbulent d issipation (Workm an &| 
Armitage 2008) Vaisala et al.||2014). Addressing questions 


related to the transport properties of the emerging turbulent 
flow,|Nauman & Blackman] (|2015|) have recently revisited the 
subject ofshear-rate dependence in MRI turbulence. The au¬ 
thors confirm the scaling of the turbulent stresses with the 
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shear-to-vorticity ratio ql(2-q ), first suggested by Abramow- 
icz, Brandenburg, & Lasota (1996), as well as the relation of 


f’essah, Chan, & Psaltis](|2006b]) quoted above. 

While the mentioned approaches primarily focus on the 
momentum equation, ideally a closure model should also en¬ 
capsulate the effects of small-scale turbulence (via the induc¬ 
tion equation) onto the long-term evolution of the large-scale 
field threading the accretion disk ( Blackman|2010) >. The rel¬ 
ative importance of the turbulent Hiffusivity versus the tur¬ 
bulent viscosity h as pre viously been assessed in unstratified 
MRI simulations (|Guan & Gammie|2009[|Lesur & Longaretti| 


2009 1 Fromang & Stone 2009). It is crucial to appreci ate that 

large-scale held acts as a net flux locally flSorathia et al. 


|20T0| and thus determines the driving of turbulence vlaTmear 
MRI - in turn affecting the local value of turbulent transport 
coefficients governing the redistribution of angular momen¬ 
tum ( |Blackman||2012) l. Beyond linear MRI, the variability 
in disk coronae may be a consequence of the e mergen ce of 


large-scale azimuthal field (Blackman & Pessah 


the long-term radial evolution of the net-vertical flux (Guilet 


and 


& Qgilvie||2012 |2013| 20 14) may have important implica- 


tions for the potentiaTdrivmg of a disk wind (see disc ussion 
in |Lubow et al.|[l994| ), although |von Rekowski et al.| ( |2003[ ) 
had found that a sufficiently strong dynamo can support out¬ 
flows even in the absence of significant vertical fields. The im¬ 
portance of me an-field effects is now also witnessed in global 
simulations (|Flock et al.|2012|| Jiang et al.|2014|) that show the 
same cyclic behavior of the azimuthally-averaged magnetic 
fields that had previously been observed in local boxes only 
(see discussion in |GresseT||2010| l. Related to this, |Stepanovs| 
et al. (12014) have recently argued that the quasi-periodic char- 


acter of the disk dynamo (which they trigger by hand) may be 
reflected in episodic ejection events, which can be observed 
as knots along the jet axis. 

As with the turbulent stresses, a promising route to nurtur¬ 
ing our insight is to understand the parameter dependence of 
the various transport coefficients. Initiated by the numerical 
result of Yousef et al. (|2008|), 
de field 


plification of large- 


who found an exponential am- 
le fields in unstratified shear flow, the 


role of shear in generating correlations in the mean electro¬ 
motive force has received renewed interest (Heinemann et al. 
|2011||Mitra & Brandenburg|2012[|Hughes & Proctor|2013| l. 

The shear-rate dependence of the dynamo a effect in MRI tur¬ 

bulence was first studied by |Ziegler & Riidiger| (|20011, who 
inferred a (pseudo-) scalar a simply as the ratio of the elec¬ 
tromotive force (EMF) to the mean field - hence ignoring 
any contributions via a tensorial a effect, and/or the turbu¬ 
lent diffusivity, n (also see Brandenburg et al.119951 Ziegler & 


Rudiger ZUUUp. me numerical study ot Ziegler & Kudiger 
was_prompted by the theoretical work of Rudiger & Pipin 
(20001, who predicted that the negative a effect (Brandenburg 


would only exist in a narrow range of the shear param¬ 
eter. We are able to confirm this prediction, but also show 
that this “buoyant” a effect is probably not responsible for 
the cyclic patterns seen in space-time diagrams of the mean 
horizontal magnetic field. For the positive a effect, however, 
we find a linear dependence on shear, which nicely explains 
the cycle period that we infer from our simulations as a func¬ 
tion of the shear rate - albeit not the direction of propagation. 
The result nevertheless lends new support to the classical crQ 
dynamo as a natural explanation of the magnetic cycles. 

The second theme of our investigation is the NVF- 
dependence of (i) the turbulent stress and (ii) the dynamo co¬ 
efficients. In contrast to unstratified NVF simulations (that 
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Figure 1 . Top: Shear-rate dependence of the Reynolds stress (squares), 
Maxwell stress (triangles), and their sum (diamonds). The dotted lines in¬ 
dicate t he sca ling with the shear-to-vorticity ratio, q /(2 - q) first suggested by 
Abramowicz et al. 1 1996j. Bottom: shear-rate dependence of the Maxwell-to- 
Reynolds stress ratio. The dotted lines indicate scaling, (4 - q)/q, introduced 
by [Pessah, Chan, & Psaltis|j2006b) . 


find a near linear scaling, Pessah et al. 2007 

) but agreeing 

with results of global stratified simulations of 

Sorathia et al. 


transport coefficients and the net-vertical flux. If one excludes 
the strong-field regime where vigorous overturning motions 
disrupt the vertical disk structure, the dynamo coefficients 
are found to be surprisingly insensitive to the amount of net- 
vertical field. 

The third subject of our paper concerns non-locality in 
the transport coefficients. In a seminal paper, Brandenburg 
& Sokoloff ( |2002| ) first pointed out the possibility of scale- 
dependent a and rj tensors, and obtained the first results 
for stratified MRI turbulence. Even though their approach 
was partly hampered by very noisy data, their derived scale- 
dependent dynamo model provided a better fit to the simu¬ 
lations in terms of the cycle period. The aspect of scale- 
dependence - or alternatively, of the non-locality in the re¬ 
lation between the electromotive forces and the mean mag¬ 
netic field - can naturally be expressed within the test-field 
(TF) framework, where a wavenumber is typically assigned 
to the TF inhomogeneity. Varying this wavenumber, Bran- 
|denburg, Radler, & Schrinner|(]2008b|l first introduced a scale- 
dependent version of the TF method. The authors have fur¬ 
thermore obtained the shape of the convolution kernel that 
expresses the non-locality of the a and ?/ effects for the cases 
of simple prescribed flows, and for forced isotropic turbu¬ 
lence. We here apply this approach to the case of stratified 
MRI simulations and obtain, for the first time simultaneously. 
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the scale-dependence of the mean-field coefficients along with 
their vertical profile. The knowledge of these functional de¬ 
pendencies will enable new, quantitative mean-field models. 

While we here restrict ourselves to an isothermal equation 
of state, it should be remarked that allowing for redistribu¬ 
tion of heat via strong convective motions (jBodo et al.|2012| 
201 3|>, the ef ficiency of the dynamo can be drastically affected 
lGressel|2013|l, potentially explaining the enhan ced turbulent 


transport related to transient accretion outbursts (Hirose et ah 
2014[ > and variabilit y in the inner regions of protoplanetary 
disks ( Hirose||2015) >. In view of devising a comprehensive 
mean-field model of large-scale field evolution, this has the 
dramatic repercussion that such a model would be required to 
encapsulate the disk thermodynamics. 

As a final introductory remark, we note that as in previous 
papers (see, e.g., Nauman & Blacfcmanj |2015| l, we focus our 
attention to the Rayleigh stable regime, q < 2. Our simula¬ 
tions beyond this value showed exponential growth of m - 0 
modes. This is expected from linear theory (Balbus 2012} 
[Shakura & Postnov|20151 l, and the loss of the local character 

of the instability clearly precludes the use of the shearing box. 

This paper is organized as follows. Section [2] describes the 
numerical method and comput ation al set-up and introduces 
the test-field concept. In Se ction [3~i~| we discuss the shear-rate 
dependence, and in Section 13.21 we demonstrate how this af¬ 
fects the cycle period of the dynamo. Section |3T3] deals with 
the NVF-dependence of th e tur bulent stress and the closure 
coefficients, and in Section [374] we derive the spectral depen¬ 
dence of the mean-field effects. We discuss our results in Sec- 
tion[4]and draw conclusions in Section^ 

2. METHODS 

The numerical simulations presented in this paper consist of 
solving the single-fluid MHD equations in the framework of 
the isothermally stratified shearing box ( Gressel|[2010 >. Im¬ 
plementation details regarding the shearing-periodic bound¬ 
aries, and the momentum source terms in the nirv ana-iii code 
are described in detail in Gressel & Ziegler ( 2007| >. 


2.1. Numerical method 

All computations are carried out using the second-order ac¬ 
curate in space, and third-order accurate in time nirvana-iii 
code (|Ziegler||2004}), which we have furthermore extended 
with a Harten-Lax-van Leer-Discontinuity (HLLD) Riemann 
solver ( (Miyoshi & Kusano|2005jl t hat uses up-winded electro¬ 
motive forces^n^^me^^^tonej^OOS ) that have been inter¬ 
polated from cell faces to edges. The use of the more accu¬ 
rate HLLD (compared with the standard HLL solver) is ex¬ 
pected to improve the accuracy for low Mach number flows. 
In the formulation with conserved variables p, p\, and written 
in local Cartesian coordinates (x, y, z) the equations of ideal, 
isothermal MHD are: 


d,p + V-(pv) = 0, 


( 1 ) 


d,(p \) + V- [pvv + p* - BB] = -pVO - 2 pQ. (zxv), (2) 


<9,B - Vx(vxB) = 0, 


(3) 


subject to the additional constraint VB = 0, which is preserved 
to machine accuracy in our constrained-transport scheme. We 
furthermore introduce the total pressure p * = p + ^B * 2 , and 
specify the effective potential 


1 


,2 2 


®(x,z) = + -Q-Z 


(4) 


in the local frame of reference, rotating with a fixed angular 
frequency Q = Oz. Note that unlike in our previous papers, 
we here absorb the minus sign into the definition of the shear- 
parameter, q = -d In D/d In r. 

For all simulations presented here, we use the same box 
size of L x x L v x L z = H x j:H x 6 H with a fixed linear 
resolution of ~ 32 /H in all three space dimensions. This 
amounts to 32 x 96 x 192 cells in the radial, azimuthal, and 
vertical coordinate directions, respectively. The smaller box 
size compared to, for instance, the recent study of Nauman & 


|Blackman| ( |2015) > - who used a box with 4 H x 4 H x 8 H at a 
resolution of 24/// - is chosen to accommodate for a larger 
number of simulations, and a longer simulation time of typi¬ 
cally 800 orbits (with period P = 2ttQ '). This allows us to 
obtain much better statistics, which is especially relevant for 
the dynamo coefficients obtained with the test-field method. 
We have furthermore performed a reference simulation with 
twice the resolution, 64///, which we ran for ~ 270 orbits. 
The studied quantities are typically converged at the 20-40% 
level (also cf. Tab. [TJ. 

We note that we consciously neglect explicit viscous or re¬ 
sistive terms (but see discussion in Fromang et al.]|2007 for 


unstratified systems) and rely on the truncation error of our 
numerical scheme to approximate the effects of small-scale 
energy dissipation. This implies that our approach is better 
suited for addressing the dynamics of accretion disks that do 
not have extreme ratios of viscous to resistive microscopic 
dissipation coefficients (jBalbus & Henri||~20081. The depen¬ 


dence of the MRI on dimensionless numbers is an active field 


of research (e.g. 

Fromang et al. 2007 iLesur & Longaretti 

2007; Pessah & C 

han 2008 Longaretti & Lesur 2010). At low 


of vertical stratification, the MRI furthermore exhibits inter¬ 


esting dynamical behavior (Herault et al. 2011 Riols et al. 
|20l3l|20l5l >. 

Even though the numerical timestep may in practice be re¬ 
stricted by the Alfven speed high up in the disk corona (espe¬ 
cially in models with a net-vertical field), we nevertheless ap¬ 
ply an orbital advection scheme (cf. |Stone & Gardiner|2010j). 
This has the advantage that the truncation error of our Go- 
dunov scheme, which sensitively depends on the advection 
velocity, remains uniform throughout the simulation domain. 
For the interpolation of non-integer displacements, we use the 
Fourier method suggested by Johansen et al/](|2009|>. Rather 


than using a so-called Alfven limiter (see, e.g. Miller & Stone 

diffu; 


2000), we instead i nclude an artificial m ass diffusion term (de- 
scribed in detail in Gress el et al.|201l] l to circumvent undue 
time-step constraints resulting from low-density regions in the 
upper disk corona. 


2.2. Initial conditions 

We initialize the fluid velocity with the background equi¬ 
librium solution Vo = qSl x y, which - in the orbital advection 
scheme - amounts to setting the perturbed velocity identically 
to zero. Neglecting the pressure support of the weak magnetic 
field, we obtain the initial density profile by solving for hydro¬ 
static equilibrium in the vertical direction, yielding 

p(z) =p 0 exp • (5) 

This equation also serves to define our adopted convention 
for the pressure scale-height, H. We start our simulation by 
seeding the density (and hence the pressure) with a standard 
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white-noise perturbation of 1% rms amplitude. 

The initial magnetic field in all the simulations has the form 


B = 



fffc) dNVF 

— +B 0 
Po 


z + B x (x,z)x, (6) 


which is designed to obtain a transition into the turbulent state 
that is as uniform as possible. By scaling the vertical field 
with a factor ( p(z )/ po )^ 2 , we obtain/? P = const (which we set 
to a value of 1600, initially), implying that the wavelength of 
the most unstable MRI mode, 4 mri, becomes independent of 
height. Mandated by the divergence-free nature of the mag¬ 
netic field, and because B z = B z (x, z), we have to add a cor¬ 
responding radial field, B x (x,z )• In practice, this is done by 
specifying a suitable vector potential A = A y (x, z ) y. In the 
cases where we are interested in the shear-rate dependence of 
the resulting turbulence properties, the initial magnetic con¬ 
figuration is of the zero-net-flux (ZNF) type, i.e., B ^ VF = 0. 
In simulations with a net-vertical field (NVF), the amplitude 
of the component Be, = By VF is typically stated in terms of the 
midplane value, /?™ d = 2po/B~, of the plasma parameter, fi v . 


2.3. Boundary conditions 

Reflecting the local character of the box geometry, the hor¬ 
izontal boundary conditions (BCs) are of the usual sheared- 
periodic type. As in previous work, we modify the hydrody¬ 
namic fluxes (and the azimuthal component of the electromo¬ 
tive force) to retain the conservation properties of the finite- 
volume constrained transport scheme (jGressel & Ziegler 
2007). In the vertical direction, we adopt stress-free BCs, 
that is, d z v x = d : v y = 0. We compute the density values 
of adjacent grid cells in the z direction to be in hydrostatic 
equilibrium (neglecting, however, any contribution from the 
magnetic pressure). We furthermore allow outflow of gas 
through the z boundaries, while preventing in-fall of mate¬ 
rial from outside the domain. In order to counter-act the se¬ 
vere mass loss owing to the emerging disk wind (even in the 
ZNF case), we continuously rescale the mass density to re¬ 
instate the original mass, Mo, within the box. This is done in 
the same way as in earlier work ( |Gressel et al.|2012) , which 
also adopted an isothermal equation of state. We remark that, 
in situations with a strong NVF, the mass loss per orbit can 
approach 100% of the mass contained within the box. Nev¬ 
ertheless, such a replenishing of material can be thought of 
as a natural consequence of radial mass transport within a 
global disk if the system is in a quasi-stationary state on aver¬ 
age. For the magnetic field, we use the standard normal-field 
d z B = 0, B x = By = 0 (sometimes referred to as “pseudo¬ 
vacuum”) boundary condition in the vertical direction. Forc¬ 
ing the horizontal field to zero at the vertical boundaries has 
the effect of creating an additional pressure gradient, poten¬ 
tially helping the already existing outflow to become super- 
Alfvenic. 


2.4. The test-field method 

One central objective of this paper is to study how the 
dynamo closure coefficients of the mean-field dynamcQ de¬ 
pend on various input parameter such as the shear-rate and 
the strength of the applied net-vertical field. As in previous 

1 These are thought to be responsible for the evolution of the large-scale 
magnetic fields via a mean-field dynamo equation - see, for instance, |Bran-| 
|denburg, Sokoloff, & Subramanian||2012j for a topical review. 


work, we heavily rely on the now well-established test-field 
method ( |Schrinner et al.||2005| |2007) to infer closure coef- 
ficients. The se describe, for instance, the classical a effect 
( Parker) 1955 i, turbulent (“diamagnetic”) pumping, y, and en¬ 
hanced resistivity, i]\, which have long been recognized to 
affect the long-term evolution of magnetic fields in turbulent 
systems ( Krause & Radler||1980) . The TF method provides 
us with time-steady measurements of the a and // tensors as 
functions of the height in the diskj^] Obtaining these coef¬ 
ficients, for instance, as a function of the shear-rate - while 
building on existing theory of mean-field dynamos - allows us 
to make meaningful predictions about whether the character¬ 
istic butterfly diagram can indeed be described by a classical 
aD dynamo, for which a dispersion relation can be derived 
(e.g., Brandenburg & Subramanian(2005). 

The method, as it is app lied here, can be classified as 
“quasi-kinematic” ( |Brandenburg et al. 2008a) in the sense that 
it only solves for magnetic fluctuations B'Qt, t). This is done 
via an additional set of induction equations solved simultane¬ 
ously with the main simulation. The measurement procedure 
is in fact similar to the approach used by|Lesur & Longaretti| 
(2009), who impose a sinusoidal magnetic field perturbation 
and infer the correlation between the resulting electromotive 
force and the current. The difference in the TF method is 
that the perturbation is virtual in the sense that the simulation 
remains unaware of the applied field. While the additional 
equations properly take into account the actual turbulent ve¬ 
locity u'(jt,f) of the main simulation, they ignore contribu¬ 
tions in the mean electromotive force, £ = u'xB' that derive 
from a non-linear evolution of the kinematic fluctuations (see, 
however, |Rheinhardt & Brandenburg|2010). Employing sim¬ 
ulations of forced turbulence, Brandenburg et ak](j2008a) have 
found that the quasi-kinematic method remains valid for cases 
that do not contain dominant magnetic fluctuations that are al¬ 
ready present in the absence of mean fields. 

The shearing-box approximation, that we use here, elimi¬ 
nates characteristic gradients in the radial and azimuthal di¬ 
rections. We accordingly consider mean fields that depend on 
the vertical coordinate only. Using the Lagrangian velocity 
u = v + qClx y, and defining u = u + u', as well as B = B + B', 
the mean-field induction equation can be derived as 

<9 f B(z) - V x [ u(z)xB(z) + £(z) - (^Qvy)xB(z) ] = 0. (7) 

This equation, via the term £(z), still depends on the fluctuat¬ 
ing quantities u' and B'. To solve it independently in the form 
of a mean-field dynamo model, one needs to parametrize £ in 
terms of the mean quantities, u, and B, and statistical prop¬ 
erties of the fluctuating quantities alone. Following Branden¬ 
burg (2008), we apply a local formulation for the EMF 


£,(z) = tt i7 (z) Bfiz) - rjij (z) £ ]Z ! d z Bi(z ), 


( 8 ) 


where the properties of the turbulence are encapsulated in the 
a,j(z ) and r/, ; (z) tensors, relating the turbulent electromotive 
force to the local values of the mean magnetic field, B(z), and 
its derivatives. 

The approach embodied in equation formally requires 
separation of scales, such that non-local dependencies in 


space, and so-called “memory effects” (Hubbard & Branden- 
bur g|2009] > can be neglected. Accounting for such effects re¬ 


quires to replace the simple multiplicative relation in equa- 


- Further details are documented in section 3.1 of Gressel et al. 2008b 
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tion (|8) by a convolution integral , as described in detail via 
equations (fl5i 


(18 i in Section 3.4 (also cf. Brandenburg 
et al. |2008b[>. In this case, the closure coefficients formally 


become co nvolution kernels. For the ca se of isotropic forced 
turbulence, [Brandenburg et al. (2008b ) have found that the 
kernels can be approximated by a Lorentzian in Fourier space, 
implying exponential decay in real space. This was also the 
case for helical turbulence with shea r (|Mitra et al.|2009]> and 
for passive scalar diffusion (Madarassy & Brandenburg|2010) 
- which may imply that thisls a generic feature. 

A test field is thought to be an external inhomogeneity 
(caused in reality by the large-scale magnetic field), in turn 
leading to non-vanishing correlations in u'xB'. In the case of 
a horizontally-periodic, stratified box, it is common to chose 
the test fields to be functions of z, and apply averages over the 
x and v directions. Despite the non-periodic boundaries in the 
vertical direction, we chose harmonic functions with a speci¬ 
fied wavenumber k TF . In our previous studies, we were mainly 
interested in large-scale mean fields, and hence restricted our- 
selves to a single value k TF = k\ = 2n/L z . In sections 
and 


3.3 


3.1 


we again focus on a single value, namely k lt = ki. 
To study a potential depen denc e on the scale-size of the mean 
magnetic field, in section [L4j we extend this approach to a 
range of values kT F = k\ ... k' 1,2 (while keeping the NVF and 
shear rate fixed for this set of runs). 


3. RESULTS 

Our simulation results are summarized in Table |T| where 
the first two columns show the input parameters, and the re¬ 
maining columns list characteristic quantities, and represen¬ 
tative mean-field dynamo coefficients. We start by discussing 
the effect of varying the shear-parameter, q, while keeping the 
net-vertical field fixed at a low value, corresponding to a mid¬ 
plane /i p - 2 x 10 5 . Conversely, we then study the dependence 
on the NVF, when keeping the shear rate fixed at its Keple- 
rian value, q = 1.5. Comparing the results for the turbulent 
transport coefficients with the inferred mean-field coefficients 
promises to provide insights into how the two phenomena are 
potentially linked. Lastly, we study potential non-local con¬ 
tributions to the mean-field dynamo tensors, keeping both q 
and /ip fixed. 
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3.1. Dependence on the shear-rate 

One criticism of the classic “enhanced” viscosity picture 
has been that the stress associated with eddy viscosity is as¬ 
sumed to depend linearly on the shear rate (see discussion 
in |Pessah et al.|]2008) . This has been shown to be an in¬ 
valid assumption, since the efficiency of the turbulent angular- 
momentum transport depends on the shear rate in a more sub- 
tle way, as has also been found by many author s (Abramowicz 


et al.||1996| Ziegler & Rudiger 2001; Liljestrom et a l.||2009| 


Nauman & Blackma n|2015| l. The upper panel of Fig. [T| illus 

trates the dependence of the turbulent stresses on the shear pa¬ 
rameter q. Our results are in excellent agreement with those of 
Nauman & Blackman| (]2015|), who recently studied the ZNF 
case, and also confirm the scaling of the total stress with the 
shear-to-vorticity ratio q /(2 - q) that had first been suggested 
by Abramowicz et al.| ([1996 ) based on their numerical find¬ 
ings! The scaling of the hydrodynamic Reynolds stress is 
somewhat steeper, leading to the known (4 - q)/q scaling for 
the ratio of the magnetic-to-hydrodynamic stresses (jPessah 
et al. 2006b), which we plot in the bottom panel of Fig.[T[ Our 
results are consistent with the stratified simulations of Ziegler 


Figure 2. Turbulent mean-field coefficients based on the TF analysis (t = 
100 - 800 orbits). As a representative example, we show the curves for the 
model with q = 1.2. The labels on the vertical axes indicate the quantities 
plotted in light (<?**, q xx , rjxy ) or dark ( ctyy ,...) shades. The mean vertical 

flow velocity, u z , in panel (b), and the quasi-linear estimate for the isotropic 
turbulent diffusion, in panel (c), are over-plotted as dashed lines. 


& Rudiger (2001) and Nauman & Blackman (2015), and the 
unstratified simulations of Pess ah et al.| ( |2008 1 and |Liljestr5m| 
et al.| (2009). Quite remarkably, this “signature” can be traced 
back to the eigenmode structu re of the linear MRI problem in 
the limit of ideal MHD (Pessah et al.|2006b) . 

In order to establish whether the saturated amplitude of the 
turbulent transport coefficients depends on the presence of a 
mean-field dynamo, we have measured the coefficients that 


characterize such activity. As in p revious work (Brandenburg 
2005 2008 1 Gressel||2010 2013) , we assume a simple local 


closure between the turbulent electr omot ive force, £, and the 


mean magnetic field, B (see section 2.4 for details), and use 
the test-field method to obtain vertical profiles of various ten¬ 
sor components. 

The results obtained for the dynamo coefficients with the 
TF method are illustrated in Fig. [2]for the model with q = 1.2, 
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Table 1 

Summary of simulation results. 


<7 

£NVF 

[/?p 800 ] 

^Keyn 

1 R<t> 

[10- 3 pol 

T-Maxw 

1 R(p 

[10 -3 po] 

'T’ratio 

1 R4> 

P eye 

[2 nQT * 1 ] 

T c 

[cr 1 ] 

M rms 

[HQ] 

peak 

a yy 

[10 2 //£T] 

buoy 

a yy 

[10~ 2 tfn] 

peak 

Tfyy 

[10^ 2 tf 2 fi] 

j-mid 

•yy 

[io- 2 # 2 n] 

0.3 

0.06 

0.01 ± 0.00 

0.24 ±0.10 

36.7 

48.1 

0.147 

0.10 ±0.03 

0.106 

-0.002 

0.091 

0.019 

0.6 

0.06 

0.07 ± 0.03 

1.13 ±0.37 

16.2 

25.2 

0.091 

0.26 ± 0.08 

0.386 

-0.011 

0.581 

0.090 

0.9 

0.06 

0.23 ± 0.09 

2.12 ±0.70 

9.3 

12.6 

0.067 

0.41 ±0.12 

0.673 

-0.036 

0.972 

0.155 

1.2 

0.06 

0.76 ± 0.26 

4.47 ± 1.35 

5.9 

9.5 

0.046 

0.68 ±0.15 

0.936 

-0.059 

1.478 

0.314 

1.5 

0.00 

1.49 ±0.46 

5.88 ± 1.63 

3.9 

7.5 

0.039 

0.86 ±0.14 

1.048 

-0.014 

1.530 

0.556 

1.8 

0.06 

4.16 ± 1.02 

11.75 ±2.88 

2.8 

5.9 

0.027 

1.31 ±0.06 

1.200 

- 

1.513 

1.232 

1.5* 

0.00 

1.14 ±0.37 

4.99 ± 1.70 

4.4 

6.8 

0.016 

0.84 ±0.01 

0.867 

-0.015 

1.233 

0.337 

1.5 

0.01 

1.49 ±0.44 

5.88 ± 1.55 

3.9 

7.3 

0.039 

0.69 ± 0.03 

1.031 

-0.017 

1.541 

0.533 

1.5 

0.02 

1.38 ±0.39 

5.48 ± 1.35 

4.0 

7.1 

0.038 

0.67 ± 0.03 

1.107 

-0.020 

1.587 

0.494 

1.5 

0.04 

1.58 ±0.55 

6.26 ±2.01 

4.0 

7.2 

0.036 

0.72 ± 0.02 

1.023 

-0.003 

1.544 

0.586 

1.5 

0.08 

1.70 ±0.53 

6.75 ± 1.99 

4.0 

7.1 

0.035 

0.76 ± 0.02 

1.100 

-0.014 

1.614 

0.603 

1.5 

0.16 

2.20 ± 0.57 

8.81 ±2.00 

4.0 

6.4 

0.030 

0.90 ± 0.03 

1.149 

- 

1.833 

0.703 

1.5 

0.32 

3.57 ± 0.96 

12.78 ±2.73 

3.6 

5.9 

0.018 

1.20 ±0.01 

0.950 

- 

2.475 

0.682 

1.5 

0.64 

5.59 ± 3.00 

15.06 ±6.17 

2.7 

- 

0.012 

1.32 ±0.03 

0.655 

- 

4.425 

0.436 

1.5 

1.28 

5.32 ± 9.49 

12.74 ± 11.33 

2.4 

- 

0.019 

1.11 ±0.04 

0.604 

- 

4.231 

0.550 


Notes: Here, the shear rate is defined as q = -d In £2/d In r, that is, q = 1.5 for Keplerian rotation. The net-vertical field, B z , is given in units that are multiples 
of the field strength resulting in a midplane plasma parameter, /?™ lcl = 800. The Reynolds and Maxwell stresses are computed as correlations of fluctuating 
quantities. The cycle period, P C yc, of the dynamo butterfly diagram is obtained as described in the text. The correlation time, r c , refers to the classic estimate 
of the turbulent diffusivity (cf. Fig. 151. In the chosen units, a mls is equivalent of a turbulent Mach number. Dynamo coefficients labeled 'peak’ ('buoy’) 
correspond to open squares (triangles)in Fig. [5] The same holds for the turbulent diffusivity, which is plotted in Fig. [4] accordingly. The run marked with the 
asterisk (*) is a double-resolution reference run with 64/H grids, indicating that results are generally converged at the"20— 40% level. 


and B z = 0.06. The first panel, (a), shows the diagonal el¬ 
ements of the dynamo a tensor. For homogeneous isotropic 
turbulence, this a effect can be shown to be propo rtional to 
the negative kinetic helicity (Krause & Radler||1980). For ki- 
netically driven, stratified, rotating turbulence, analytic theory 
predicts that a will emerge naturally from the inhomogeneity 
introduced by the stratification combined with the anisotropy 
caused by the rotation. The magnitude of the effect is pre¬ 
dicted to be proportional to the gradient in the mean density 
combined with some power of the gradient in the turbulent 
velocity (see the discussion in Brandenburg et al.|2013J>. 

Because of the dominant shear, a xx is expected to play a 
minor role, since azimuthal field can efficiently be generated 
from radial field by the background shear flow - the so-called 
ail dynamo regime. Then, the a yy element plays an important 
role in producing radial fields from azimuthal fields via the 
vertical variation of £> y = a yy B y . Note that a yy shows opposite 
slopes near the disk midplane and in the disk corona. 


Based on the known behavior of ail dynamos, Branden¬ 
burg ( |l998| l has argued that a negative slope near the mid¬ 


plane is required to explain the upward migration observed in 
the characteristic “butterfly” diagram in stratified MRI simu¬ 
lations (also see discussion in|Gressel|2010|l. The results pre¬ 
sented here challenge this hypothesis, since we find the neg¬ 
ative slope to vanish for q ^ 1.5, whereas t he dy namo pattern 
remains unaffected (see discussion in Sect. |3.2| >. This finding 
potentially implies that another effect may be responsible for 
the observed pattern speed (and direction of propagation). 

A peculiarity of the MRI dynamo coefficients that had been 
noted before (jBrandenburg||2008) can be seen in the next panel 
of Fig. panel (b) shows the off-diagonal elements of the 
a tensor, which are found to be symmetric (i.e., a xy — a yx ). 
In isotropic kinematically-forced turbulence, one typically ex¬ 
pects (e.g., |Krause & Radler|1980) that the dominant contri¬ 
bution to the off-diagonal elements are given by so-called dia¬ 
magnetic pumping, y = -V;/-| /2, reflected in auh-symmetric 


tensor elements, a yx = -a xy . Such an effect implies an 
expulsion of magnetic flux from strongly turbulent regions 
into more quiescent regions, and has for instance been ob¬ 


served in kinetically forced interstellar turbulence (Gressel 
|et al.||2008a] l. The fact that both coefficients show the same 
sign implies that this effect is significantly anisotropic in MRI 
turbulence, that is to say that B x and B y are affected differ¬ 
ently. 

The diagonal elements of the diffusivity tensor are shown 


in panel (c) of Fig. |2| Complementary to the result of Lesur 
& Longaretti ((2009T, who find the diffusion tensor to be 


anisotropic when comparing the radial diffusion of the az¬ 
imuthal and vertical field components, the coefficients mea¬ 
sured from our simulations approximately satisfy rj xx - q yy 
- these coefficients however describe vertical diffusion of the 
azimuthal and radial field components, respectively. Isotropy 
is not necessarily expected for a shear flow with a dominant 
azimuthal field component (Vainshtein et al.j [1993 1, and the 
deviation is indeed most pronounced in the more st rongly 
magnetized disk corona (z <; 2 H). As previously found ( Lesur| 
|& Longaretti|2009~||Gres sel 2010), the q tensor, moreover, has 
significant off-diagonal contributions, shown in panel (d). The 
r] yx term has the potential to couple radial and azimuthal fields 
(via their vertical gradients), and hence may be responsible 
for pattern propaga tion and/or cyclic behavior (Gressel 2 010[ 
|Herault et al.||20lT] l. The three dominant coefficients of the 

diffusivity tensor (i.e., q xx , rj yy and q xy ) closely follow the the¬ 
oretical expectation 


1 9 

V “ 3 t c M rms (9) 

within the disk body (z ^ 2 //, see panel (c) in Fig. [2]), where 
we have assumed a constant correlation time, r c , as well as 
isotropic velocity perturbations. In the low -[5 disk corona, 
magnetic fluctuations are likely dominant, perhaps explain¬ 
ing the deviation from the kinematic result. We will later use 
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Figure 3. Vertical profiles (main plot), and shear-dependence (inset) of the 
a yy (z) coefficient. Amplitudes for the kinematic dynamo (squares), and the 
buoyancy-driven dynamo (triangles) have been inferred from peak values as 
indicated. Heuristic approximations for the ^-dependence are also shown 
(dashed lines). When approaching the Rayleigh limit (g—>2), field-line buoy¬ 
ancy appears to be more and more dominated by hydrodynamic effects. 


equation (|9| to obtain an independent estimate of the turbu¬ 
lent coherence time r c as a function of shear rate (see Fig. [5] 
below). 

In Fig. [3] we compare the vertical profiles of the a yy co¬ 
efficient for the six simulation runs with varying shear rate J^] 
For lack of a better alternative, we resort to describing the 
strength of the effect simply by peak amplitudes (squares). 
The dependence of these amplitudes is plotted in the inset of 
Fig.[3](again using open squares), where we see that the dom¬ 
inant positive a effect scafes approximately linearly with the 
shear rate. The shear-rate dependence is markedly different 
for the weak negative a effect near the disk midplane (trian¬ 
gles). This effect is thought to originate from correiated mag¬ 
netic fluctuations ( )Brandenburg| 1998 1 that are associated with 
the buoyant rise of reg ions with strong azimuth al field. Based 
on quasi-linear theory, Rudiger & Pipin ( j2000 ) had predicted 
that this effect should only be present for a finite range of the 
shear parameter - given by the condition 0.75 < q < 1 in 
their analysis; cf. their equations (14) and (20). While the 
precise parameter range is different from this analysis by a 
small factor, we indeed confirm this predicted behavior. We 
empirically describe the q dependence to be proportional to 
<gr 4 (1.5 — q), and we speculate that the discrepancy may be re¬ 
lated to the assumed spectrum of the magnetic fluctuations in 
their derivation; cf. |Riidiger & Pipin| (j2000|l, equations (17) 
and (18). 

The result illustrated in F ig, [^qualitatively agrees with that 
of Ziegler & Riidiger[ ( |2001] >, who also found a to vanish at 
a finite value of q, and whose values scale similarly to ours - 
albeit their effect is stronger by about a factor of four (cf. their 
table 1). In general terms, the existence of a finite interval for 
this effect derives from the individual contributions due to the 


shear and Coriolis term entering with opposite signs (Rudiger 
|& Pipin|2000) l. We note that this is however qualitatively dif¬ 
ferent for the positive a effect further up in the disk (which 


may be more akin to the mechanism derived in |Riidiger & 


3 To facilitate a direct compariso n, we have filtered the profiles by expand¬ 

i ng them into l.cgcndrc polynomials up to order 12, and set the even (od(T) 
components to zero for a (q) to cnlorce thc analytically expected symmetry. 



-0.5 0.0 0.5 1.0 1.5 2.0 2.5 5.0 

height l [H] 


Figure 4. Same as Fig. [3] but for the diffusivity coefficient // vv (V). Peak 
values in the disk coronafsquares) roughly scale linear with q, while the 
midplane values show a deviant scaling with q ~, instead. Qualitatively similar 
trends are found in the other three tensor coefficients. 



0.5 1.0 1.5 

shear rate q 


Figure 5. Shear-rate dependence of the correlation time r c (squares) derived 
from the diffusivity via assuming r\ yy - r c /3 iq ms . The solid line shows an at¬ 
tempt to independently obtain r c from the auto-correlation function of kinetic 
helicity. For comparison, we also show the scaling of u rms itself (triangles), 
which scales approximately with r/ 3/ -. explaining the quadratic scaling of 
Vyy (z = 0) seen in Fig.|4]above. 


|Kitchatinov|1993[ ). 

The vertical profiles and amplitudes of the turbulent diffu¬ 
sivity i]yy(z) are shown in Fig. [4] We again distinguish be¬ 
tween peak values in the upper disk layers (squares) and val¬ 
ues near the midplane (triangles). The former effect scales 
roughly with the shear rate, q, whereas the latter scales with 
q 2 , instead. This scaling can be understood in terms of the 
scaling of r c , and w rms , respectively. We plot both relations in 
Fig. [5] and we find r c ~ l/q, which is in goo d agreement with 
the findings of Nauman & Blackman ( [2Q15} >, who derived this 
from velocity auto-correlation functions. We moreover find, 
Mrm S ~ q -F 2 , and together this is equivalent to the q(z-0) ~ q 1 
dependence shown in Fig. [4] Since we obtain r c from the as¬ 
sumed relation (|9|, the two are not independent, of course. 
We hence independently derive a scaling of r c from the auto¬ 
correlation function of the kinetic helicity /?kin = u'- (Vxu'), 
which we use as a rough proxy for velocity correlations. This 
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Figure 6. Dynamo cycle period, P cvc a s a function of the shear parameter. 
The periods measured via equatio n (13) are shown as diamonds, and match 
well the prediction from equation d 10b, with peak values for taken from 
Table[j] and k z = 0.33. The fitted relation ( 11 } is shown as a dotted line. 

is plotted as a solid line in Fig. [5] and exhibits a similar scaling 
behavior. 


3.2. Cycle Period as a Function of Shear Rate 

In the shear-dominated limit, \a\k z (q£lf l <s 1, the cycle 
frequency of the akl dynamo can be approximated as: 

1/2 


^CVC “ 


2 a yy 


( 10 ) 


see, e.g., eqn. (6.40) in |Brandenburg & Subramanian|((2005jl. 
Expressed in terms of the cycle period, P cyc = 2jrJco cyc , and 
including the scaling relation for a from Fig. [3] we predict the 
cycle period of the butterfly diagram to scale a^] 


2 n 

p ~ _ 

cyc £1 


n q 


aoq k-H 
Jm ~2:n 


-1/2 


11.0 2 n 
q O 


~pr - (ID 


where, in the last step, we have used a constant o-q - 
0.008 HQ., corresponding to the dashed line labeled q\ in 
the inset of Fig. [3] and have arbitrarily assumed a wavenum¬ 
ber k z - 0.33 in units of 2 n/H. 

It is not straightforward to extract a meaningful wavenum¬ 
ber from the vertical profiles of B x and If , but since expres¬ 
sion ([Toll a l so appears in the real part of the eigenvalue. 


<Re T ± ^ — 7/yy k 2 + 


q a yy 


1/2 


( 12 ) 


cf. eqn. (6.39) in Brandenburg & Subramanian (2005), there 
is an alternative way to obtain k z . Using our valuesTor q yy (see 
Tab. [TJ, and assuming that the dynamo is one-dimensional 
( k 2 = k 2 ), and marginally excited CRe A = 0), we can use 
this to determine an independent estimate for the wavenum¬ 
ber, which we find to be k z - 0.5, using peak values for q yy . 
Consistent with the scalings we find for a and i], and the way 
these quantities appear in equation ( | 1 2[ >, this value is largely 
independent of the shear rate, q. 

To reliably extract the cycle period from the space-time but¬ 
terfly diagram observed in the various simulation runs, we ap¬ 
ply the following procedure. We take B y (z,t) profiles, that 
have been obtained by averaging in the x and y directions, 

4 In view of the eleven- year solar cycle, the unique numerical constant of 
eleven ‘years' in equation jilt is a rather ironic coincidence. 


and normalize them with the maximum at each given time t. 
This is to artificially enhance the contrast of the cyclic pattern. 
We then compute the autocorrelation function of I fit, z = z,) 
at each height with |z,| > 2.25 H. This exploits the fact that 
the cycle is seen most pronounced in the upper disk layers. 
Since B y (t) is periodic in time, so is the ACF, B v (t), and we 
apply a least-square fit to model the ACF as a mixture of an 
oscillation and a purely stochastic part, that is, 

B y (j) = [(1 - a) + a cos(27rcu cyc r)J e~ T/Tc , (13) 


with the coefficient a e [0,1] being the relative amplitude 
of the oscillatory part. A function of this kind was used in 


Gressel et al. (2011) for a different purpose. We find that this 
provides a useful method to extract the quasi-periodic cycle 
frequency u ) cyc , and we obtain cycle periods of P cyc = 48.1, 
25.2, 12.6, 9.5, 7.5 and 5.9 orbits for the models with q = 
0.3, 0.6, 0.9, 1.2, 1.5, and 1.8, respectively. These values 
are plotted together with the prediction ( fT0| > and the derived 
scaling (fTT|) in Fig. [6] as a function of the shear parameter. 
Given the crudeness of the approach, the agreement is striking 
and encourages further investigations into this direction. For 
instance, equation ( p~0] > assumes a constant a coefficient and 
neglects any spatial variation. Moreover, we have taken the 
liberty to specify the wave number, k z , of the dynamo mode 
ad hoc , roughly based on the shape of the B y (z) profile. 

Another reason for which our findings are surprising arises 
from the fact that we have used the values of the dominant 
positive a effect (i.e., the squares in Fig. R), and not the neg¬ 
ative effect (triangles) near the disk midplane as suggested 
by Brandenburg (1998) to explain the ‘poleward’ migration 

>□ 


of the dynamo patternrl This notion is further supported by 
our observation that the negative, buoyancy-related effect van¬ 
ishes for q 1.5, whereas the migrating dynamo pattern per¬ 
sists for all studied values of the shear parameter. If the ev¬ 
idence obtained via the TF method holds up, this means that 
the explanation for the upward propagation of the field pat¬ 
terns will have to be revised. This should however not be seen 
as a failure of the mean-field dynamo framework. After all, 
the notion that a positive (negative) a effect leads to equator- 
ward (poleward) migration of field belts in the aQ dynamo, is 
based on the assu mption of a constant isotropic a tensor (e.g., 
Brandenburg & Subramanian [2005jl. This is not the case for 
MRI-driven turbulence, where we observe a strong symmetric 
contribution in the off-diagonal elements of both the a and // 
tensors. This will mandate a thorough derivation of the mean- 
field dispersion relation, considering the mutual strength of 
the closure coefficients as well as their spatial variation. 

Figure [7] shows the prediction of the pattern speed from the 
most simple aQ model with a constant, isotropic, a effect. 
Because of the “wrong” sign of our a, the propagation di¬ 
rection has been chosen arbitrarily. Note that the strict rule 
about downward (upward) propagation for positive (negative) 
a was in fact derived from a homogeneous model. In an inho¬ 
mogeneous situation such as ours, the dynamo can be affected 
by other factors - see, for instance, Brandenburg et al.](|2009|), 
who found a reversal of the direction for wave propagation for 
different vertical boundary conditions in a stratified dynamo 
model. On top of the effect that is caused by the advective 


5 Using the negative a effect did not provide a good fit with the cycle data, 
mostly because the amplitude is too low to match the order of magnitude with 
any reasonable choice of k z . Note, however, that Brandenburg & Sokoloff 
1 2002 1 found the cycle period to be affected by non-locality in a. 
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time t [2n (qCl) -1 ] 350 


time t [277 (qQ) '] 350 


Figure 7. Azimuthal-field butterfly diagram for the models with q = 0.6, 0.9, 
1.2, 1.5, and 1.8 (top to bottom). Note that the time axis has been rescaled 
with q to allow for a more direct comparison. The black lines illustrate the 
the cycle period and pattern speed, c = (o cyc kZ 1 (albeit with an arbitrarily 
chosen propagation direction), of the aCl dynamo wave, and also include the 
advective contribution u z , which becomes apparent for |z| S 2 H. 


contribution, u z (responsible for the curvature of the charac¬ 
teristics shown in Fig. |7J, there is a slight tendency towards 
faster propagation at larger z. This illustrates that a more de¬ 
tailed dynamo model is required to explain the field evolution 
in more detail (also cf. figure 6 in Gressel|2010| . 


3.3. Net-vertical field dependence 

Guided by the endeavor of developing a mean-field frame¬ 
work relevant to turbulent transport of angular momentum in 
accretion disks, we now turn our attention to the influence of 
varying levels of net-vertical magnetic flux. 

It has been understood very early that the amount of net- 
vertical field threading the disk has a direct influence on the 
accessible linear MRI modes even in the state of vigorously 
developed turbulence (Hawley et al .|l995j >. The wavelength, 
Tmri, of the fastest-growing linear models directly related to 
the Alfven speed corresponding to the vertical field compo¬ 


wavelenqth X k , RI /Ax 

o.i i.o y 10.0 



Figure 8. Net-vertical field dependence of the turbulent Reynolds stress 
(squares), Maxwell stress (triangles), and their sum (diamonds) for fixed 
q = 1.5. The horizontal dotted lines indicate ZNF values, while the slanted 
lines follow a ~ dependence. For the strongest field case, the disk is ef¬ 
fectively disrupted by the vigorous MRI turbulence. The stress ratio is found 
to be independent of the NVF strength. 

nent, that is 




32 


15y8 p 


mid 


H , 


for q=~. (14) 


For the case of unstratified MRI, Pessah et al.| ( j2007] l have 
found a linear scaling of the turbulent stresses with the most 
unstable wavelengtlff /I mri, and accordingly the amplitude of 
the net-vertical field/}.. The authors, however, point out that 
this scaling is likely to change when accounting for the buoy¬ 
ant loss of magnetic flux, which is likely unavoidable in the 
stratified case. In a previous paper that included a magneti¬ 
cally decoupled midplane in apoorly ionized, stratified pro¬ 
toplanetary disk (Gressel et al.|2012[l, we found the Maxwell 


stress to scale oc but this precluded the effect of turbu¬ 
lence near the midplane. 

In Figure[8] we plot the dependence of the turbulent stresses 
on the net-vertical field strength. The influence of the NVF 
becomes noticeable at /i|’" d - 10 5 , and the additional contri- 
butior| 7 |to the turbulent stresses roughly scales ~ /j p 1 4 . The 


break (Pessah et al.|2007|i seen in Fig. |8] furthermore appears 
to coincide with the ability to resolve theMRI in the disk body 
(see the upper axis of Fig. [8] where we show the wavelength, 
d MRI , °f the fastest-growing linear mode in units of the ra¬ 
dial grid spacing. Ax). We remark that, because /3 p < /}™ d , 
this resolution indicator provides a lower bound in the strati¬ 
fied disk. As such, the transition point between ZNF and NVF 
may shift further towards the left if going to higher resolution, 
but attacking this question is beyond the scope of the current 
paper (however see Sect. 3.3 in Sorathia et al .|2010| for a po¬ 
tential explanation). The slope of ~ /3 p ' 4 seen i 


1/4 in Figgis 


6 Note that the simulations considered in Pessah et al. .2007, correspond to 
those carried out by Sano et al. 1 2004), who reported a scaling proportional to 
^MRr This steeper scaling is derived when the simulations with vp^l{L z Q) = 
6.25 x 10 -6 and 1 x 10 -4 are excluded in the fit (see Figure 2 in)Pessah et al.|| 
1 2007). Including these results leads to the linear scaling previously reported 
by [Hawley et al.j{1995||1996). 

' This interpretation tacitly assumes that the NVF case is bounded from 
below by a ZNF self-sustaining state, whose convergence with numerical res¬ 
olution is however still under deba te - even in stratified shea ring-boxes that 
support mean-field dynamos pavis et al.|2010[|Bodo et ahj2014 
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siderably less steep than for unstratified simulations, and we 
attribute this potentially to the importance of buoyant flux loss 
and the modified eigenmodes f or stratified case ([Latter et al. 


2010[ McNally & Pessah|2014j l. The observed flux-stress re¬ 


lation is qualitatively simila r to the one seen for in termediate 
disk layersn in Figure 2 of Sorathia et al. (2010 1 , who per¬ 
formed global ZNF simulations and analyzed whether small 
local patches of these were compatible with the flux-stress re¬ 
lation obtained in NVF box simulations. 

In our model with the strongest net flux, we observe notice¬ 
able disruption of the disk. This is in agreement with simi¬ 
lar simulations by Bai & Stone ( 2013) >, who found magnet¬ 
ically dominated disks (i.e., < 1 throughout the domain) 

for models with initial < 10 3 . In our models, the strong 
turbulence goes along with severe mass loss through the up¬ 
per boundaries (which we artificially replenish by rescaling 
the mass density in the bulk of the disk). The mass loss rate 
is on the order of Mq (2n) 1 Q, where Mq is the mass initially 
contained in the simulation. Because of the strong outflow, 
simulations with taller vertical domains than +3 H, as used 
here, should be performed (cf. Fromang e t al,|2013) . 

Keeping in mind these limitations of our current model, 
we nevertheless proceed with our investigation about how the 
NVF potentially affects the mean-field dynamo. In Fig. [9j 
we show a representative example of dynamo profiles for a 
moderately strong net-vertical flux case. In comparison with 
Fig. [2] for the ZNF case, we observe a significantly enhanced 
a xx effect, whereas a yy appears to be largely unaffected. Be¬ 
cause of the dominant role of the q£l B x term in the mean-field 
induction equation, even a xx - 0.1 still provides only a rel¬ 
atively weak contribution to the overall a 2 £l dynamo. The 
coefficients a xy , and a yx are less similar, and the latter also 
shows stronger fluctuations, as indicated by the shaded area. 
A further difference in the presence of a net-vertical field, is 
the disk wind, u z that already rises steeply near the disk mid¬ 
plane. 

Unlike in Fig. [2] the turbulent diffusivities rj xx , and j/ yv 
now show a significant deviation from the kinematic result, 
prompting at the importance of magnetic fluctuations even 
for |z| ^ H. Alternatively, the correlation time, r, could be¬ 
come dependent on the vertical position in the disk. Simi¬ 
larly to a xx , the off-diagonal component r\ xy is enhanced in 
the presence of a strong net-vertical field, even though, again 
the much smaller component rj yx is the more interesting one 


in its relation to the dominant shear (cf. appendix B in Bran- 
denburg||2005|). Despite the pronounced differences that 


have just illustrated, the actual scaling behavior of the rele¬ 
vant dynamo coefficients with the NVF is found to be much 
less dramatic. 


The dependence of a yy (z) on /j"" d is illustrated in Fig 


10 


where we plot their vertical profiles and infer the associatec 
amplitudes. This is again done separately for the region near 
the midplane (triangles), and the upper disk layers (squares). 
The impact on the weak negative a effect is the most severe, 
and the effect vanishes for any appreciable amount of net- 
vertical flux (cf. Fig. [81. Compared with the shear-rate de¬ 


pendence shown in Fig|3j the variation of a yy (z) in the upper 
disk with the NVF is minor, and only for very strong fields 
we observe a mild quenching. We speculate that this effect 
is related to the stronger outflow in this case, removing mag¬ 
netic fluctuations before they can develop significant correla- 


8 That is, at z = ±2 H, depicted by green circles and blue dots, respectively. 






z [H] 


Figure 9. Same as Fig. [5] but for the run with q = 1.5, and B z = 0.32, 
corresponding to a midplane /J p - 7800. 


tions with the fluctuating velocity. A quantitative assessment 
of this amounts to comparing the rotation timescale, ~ £1 ' 1 , 
with the advective removal timescale, ~ // it 1 . We remark 
that the observed trend is not compatible with what would be 
expected from an increasing effect of the advective helicity 
flux (Brandenburg et al. 2009), that would likely be more pro¬ 
nounced in the case of a stronger NVF (and hence a stronger 
outflow). By re moving small-scale curr ent helicity from the 
system (also see|Oishi & Mac Low|2011]l, such a helicity flux 
would allow the dynamo to “un-quench’ itself. It is however 
not straightforward to predict how such an effect would be 
reflected in the inferred TF coefficients that we present in Fig¬ 
ure [To] 

The influence on the turbulent diffusivity, when applying a 
NVF is shown in Fig. 


0 ), plotted with triang 


11 where the midplane value, rjyyiz = 
es, shows a similar quenching charac¬ 


teristic as seen for c%(z). Unlike the dynamo effect, the dif¬ 
fusivity displays an anti-quenching behavior in the upper disk 
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Figure 10. Vertical profiles (main plot), and net-vertical field dependence 
(inset) of the filtered a yy (z) dynamo coefficient. As in the corresponding 
Fig. [3] amplitudes for the kinematic dynamo (squares), and the buoyancy- 
driven dynamo (triangles) have been inferred from peak values as indicated. 


layers (squares). Given the dominant role of magnetic tension 
forces (y8 p ^ 1) in the disk corona this is not necessarily sur¬ 
prising, and it fits together with the observation that the T) yy (z) 
profiles deviate from the kinematic estimate in this region. 

Based on the simple estimate ( p~0] >, which predicts the dy¬ 
namo cycle frequency to dependon the square-root of the 
a effect, and looking at Fig. [TO] we would not expect a strong 
variation of the observed cycle period with varying NVF. Un¬ 
fortunately, the unambiguous determination of the cycles is 
hindered by the vigorous fluctuations in the simulations with 
significant net-vertical flux. A similar trend was reported 
by |Bai & Stone| (|2013jl who found more irregular cycles for 
strong NVF, and state that the dynamo is “fully quenched” 
for strong background fields (i.e., /3™ d = 10 2 , where indeed 
no cycles are observed). With the rather moderate quench¬ 
ing seen in Fig. [TO] we would however not interpret this as 
the absence of cyclic behavior (and consequently a mean-field 
dynamo) but merely as the cycles being dominated by other 
effects, such as linear MRI (see discussion in Gressel 


and parasitic modes (Goodman & Xu 1994] Latter et al. 


20101 


2009 


|Pessah & Goo dman 2009, |Pessah||20lU] C which are more- 

over expected to develop on comparatively larger scales in 
the strong NVF regime. In contrast to the findings presented 
here. Figures 10 and 11 in |Gressel et al.| ( |2012] > show a pro¬ 
nounced dependence of the cycle period on the applied net 
flux. This was however for a partially ionized protoplane¬ 
tary disk model, where the ability of the MRI to develop was 
severely limited by significant Ohmic dissipation in the disk 
body. As mentioned before, the flux removal timescale likely 
plays a role in the presence of the strong wind, and this is of 
course not reflected in equation ([TO]). 


3.4. Scale-dependent closure coefficients 

The functional dependence proposed in equation ([8]) as¬ 
sumes that the EMF at a given point in time and space can be 
modeled in terms of the mean magnetic field, and its deriva¬ 
tives, evaluated at the same point in time and space. This 
assumption can be relaxed so that the EMF is no longer local 
in space, according to 

£i(z) = J' &ij(z, 0 Bj(z~0 - f)ij(z, 0 e Jz i d z B t (z-0 dZ. (15) 



height z [H] 


Figure 11. Same as Fig. 1 10| but for the diffusivity coefficient rj yy (z). Peak 
values in the disk corona (squares) are enhanced at high NVF strength, while 
the midplane values appear to be quenched instead. Qualitatively similar 
scalings are found for the other tensor coefficients. 


Here, the tensor components &ij(0 and ffrjiO play the role of 
convolution kernels that may additionally depend on the ver¬ 
tical coordinate z - reflecting the changing conditions with 
height in our stratified model. The spatial scale £ over which 
they decay provides a measure of how the EMF at a given lo¬ 
cation depends on its surroundings. Note that equation ([15 i 
recovers d8} if ay(z, O - a,j(z) 6(0 and %(z, O = dijiO <fe), 
where <)(()ts Dirac’s delta function. In the following, we will 
omit the z dependence in our notation for clarity. Because 
the EMF is a convolution in real space, its Fourier transform, 
denoted with a tilde, is simply 


— &ij(k z ) Bj(k z ) i]ij(k z ) i k z Sj Z i B;(k z ), 


( 16 ) 


which provides the context for considering scale-dependent 
closure coefficients and illustrates the close relation between 
scale dependence and non-locality (with conjugate variables 
£ k-). 


jBrandenburg & Sokoloff] (]2002|) first mentioned the possi 
ility of a scale-dependent, or non-local, a effect in the con 


bi 

text of stratified MRI (cf. their Figures 10-12). Even though 
their approach was affected by poor statisticsfjthe overall as¬ 
sumption appeared very justifiable. A more systematic ap¬ 
proach to a non-local formulation for the EMF via scale- 
dependent aij(k z ) and ijij(k z ) coefficients can be naturally aided 
by employing scale-dependent test-fields. When the test fields 
appear as an inhomogeneity with a specific vertical wavenum¬ 
ber kl v (Brandenburg 2008) this allows us to characterize the 
response of the system as a function of this imposed scale. 
By surveying a number of scales we can measure ajj(k TF ) and 

*7#D- r _, 

Based on a series of harmonic test fields,[Brandenburg et al. 
(j2008b) first introduced a scale-dependent TF method, and 
determined the shape of the convolution kernel for isotropi¬ 
cally forced homogeneous turbulence. They showed that the 
Fourier transforms of the convolution kernels can be charac- 


9 Lacking the additional bits of information that the TF method so ele¬ 
gantly provides, their inverse problem was overconstrained and they had to 
resort to making assumptions about the tensorial structure of the T) coefficient. 






















































12 



height z [H] 


Figure 12. Vertical profiles (main plot), and wavenumber dependence (inset) 
of the a yy (z) dynamo coefficient. As in the corresponding Fig. [3| amplitudes 
have been obtained for the kinematic dynamo (squares), and me buoyancy- 
driven dynamo (triangles). While the former scale similar to a Lorentzian 
(dashed) with k c = 3.7, no clear scaling is seen for the latter. 


terized by Lorentzians, 


a(k z ) = 


Q'o 


i + (k z /k c y- 


n(kz) 


no 


i + (k z /k c y- 


(17) 


This implies that the convolution kernels in real space are de¬ 
caying exponentials. 


a(0 = y exp (~k c \£\), f/OT) 


no 


exp(-£ c |£|), (18) 


where the characteristic wavenumber, k c , can in principle take 
individual values for a and rj. In fact, Brandenburg et al. 
( |2008b[ ) found the characteristic scale over which the diffu- 
sivity 7) falls off to be about a factor of two smaller than for a , 
that is, k ( y ] =* 2 k^\ 

In what follows, we extend this approach and apply it to 
stratified MRI turbulence. While it is in principle possible to 
obtain all wavenumber components simultaneously, we have 
chosen to run a separate model for each wavenumber k TF . Un- 
like for the h omogeneous turbulence studied by Brandenburg 


et al. 


(2008b), our coefficients a(z, kJ F ), and q(z, kJ F ) are SF 
multaneously functions of the position z, and the wavenumber 
k VF . The latter, however, only enters via the prescribed test- 
field inhomogeneity (as indicated by the label ‘TF’) in the 
evolution equation for the passively evolved magnetic fluctu¬ 
ations. While this may appear somewhat confusing, we point 
out that the variables z, and k TF are not complementary in 
the usual sense, that is, when applying a Fourier transform to 
some function f(z) to obtain f{k z ). 

Accordingly, in Fig. fl2| the ordinates in the main plot and 
inset are just as independent as in all the previous plots, and 
it is meaningful to look at the k z dependence at various posi¬ 
tions z- The picture for the weak negative a effect (triangles) 
is uncertain, and no clear trend is seen with wavenumber. The 
stronger, positive a effect is flat for low wavenumbers, and 
its spectral shape is compatible with a Lorentzian of width 
!» - 
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k " = 3.7 in units of 2i t/H (note that the axis in Fig. 
in units of 2 n/L z ). If we look at Fig. [13] a similar scaling" is 
found for the turbulent diffusivity, q yy (z) - both for the disk 
corona (squares) for which we find k^ = 3.3, and for the 



Figure 13. Same as Fig. |12| but for the diffusivity coefficient rj yy (z). Peak 
values in the disk corona (squares) as well as the midplane values (triangles) 
are both compatible with Lorentzian shapes (dashed lines), 
midplane values (triangles) for which we find k F,) = 4.8. As 
we have mentioned already, [Brandenburg et alT ( [2008b i had 
found a wider kernel (by a factor of two) for a, compared to 
rj in their simulations. Here we find the opposite trend with 


k^ being roughly 15% larger than k P when comparing the 
results obtained at z — 2.5 H. However, the value of k { c n> in¬ 
creases by about 50% near the midplane, likely reflecting the 
shorter correlation length at higher density. Compared with 
the wavenumber of the large-s cale dynamo mode, that we es¬ 
timated as k- - 0.33 in Section 3.2| a value of k c - 3.3 implies 
that our simulations have a scale-separation ratio of about ten. 
While this is sufficiently large to warrant a mean-field descrip¬ 
tion in general, it may not be comfortably so - hence making 
it worthwhile to consider non-local contributions to the EMF 
in addition. 


4. SUMMARY AND DISCUSSION 

We continue our endeavor to build a unified framework 
for angular momentum transport in accretion disks, on one 
hand, and the redistribution and amplification (via a dy nam o) 
of large-scale magnetic fields, on the other hand ( Blackman] 
). In the present paper, we have studied various param¬ 


eter dependencies of the turbulent transport coefficients. The 
main focus was on the shear rate, q, and the net-vertical field 
strength (quantified by the midplane plasma parameter /3 p ). 
Simultaneously, we have investigated how the closure coeffi¬ 
cients for the classical aCl dynamo depend on these parame¬ 
ter variations. The correspondence between these two aspects 
promises to increase our knowledge about the inner work¬ 
ings of accretion disk dynamos and ultimately about the mass 
transport within such disks. 

Following a complementary avenue, we have studied the 
scale-dependence of dynamo closure coefficien ts, gr eatly im¬ 
proving over a previous attempt by |Brandenburg & Sokoloff[ 
(2002). Compared to their result, the TF method offers a much 
higher fidelity in determining the scale-dependence of the dy¬ 
namo coefficients, and at the same time allows us to simulta¬ 
neously obtain;;-dependent profiles. 

4.1. Shear-rate dependence and dynamo cycles 

Our results on the shear-rate dependence are congruent with 
recent work by Nauman & Blackman|(2015), who found that 
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the total (kinetic plus magnetic) turbulent stress scales with 
the shear-to-vorticity ratio, q/(2 - q), as previously suggested 
by Abramowicz et aL] (1996}. As predicted by |Pessah et al. 
(2006bI, and as also found by Nauman & Blackman (2015), 


.Reyn 


the Maxwell-to-Reynolds stress ratio scales 7’j^ axw /T i 
(4 - q)/q, and we find a value of - 4 for Keplerian shear. 

As in previous work (|Gressel||2010| (2013}, we apply the 
test-field method to obtain measurements for the mean-field 
a effect, and eddy diffusivity, q, which we find to significantly 
depend on the background shear-rate. The turbulent diffu¬ 
sivity is found to scale like q ~ q * 2 * * near the disk midplane, 
whereas a scaling q ~ q is found in the upper disk layers. The 
quadratic scaling can be explained by the dependence of the 
turbulent correlation time r ~ l/q (also cf. Nauman & Black- 

2015}, and the scaling of the turbulent velocity dispersion 

' 3/ 2 


u n 

M rms H • __ 

Motivated by the claim ( (Brandenburg et al.|j 1995| ) that the 
‘poleward’ migrating dynamo pattern requires a negative a ef¬ 
fect, we separately inspected the negative a yy that we find near 
z — 0.5 H , and the positive a yy measured in the disk corona 
(z — 2,5 H ). The former effect was predicted by Rudiger & 


Pipin (|2000J) to only exist in a narrow window of shear-rates, 
q. We indeed confirm this behavior by fitting a heuristic func¬ 
tion cfi (1.5 - q) . In contrast, the stronger, positive a effect 
was found to exist for all studied shear rates q = 0.3,..., 1.8, 
with a linear scaling a yy ~ q. A linear scaling with shear has 
been found previously in a study of unstratified turbulence 
subject to artificial forcing (Yousef et al.|2008}. In their case, 
the dynamo growth rate was found to depend linearly on the 
shear, as explained by linear theory (Heinemann et a l,|2011} . 
We remark, however, that this effect is likely distinct from 
the present case, where the dynamo effect is probably enabled 
by stratification rather than shear alone (see the discussion in 
Brandenburg et al.|2013|). 

The above interpretation tacitly assumes that the observed 
scaling is a result of the production mechanism. Since the 
measured a effect is likely in a non-linear, saturated regime, 
the scaling may equally well reflect the quenching mecha¬ 
nism. This opens the possibility to test ideas related to shear- 


driven helicity fluxes (Vishniac & Cho 

2001 [ Brandenburg 

2008; Vishniac||2009 Hubbard & Branc 

enburg |30l 1, Vish- 


niac & 


Shapovalov||2014} 

aturation. Regard!* 


, and their role in determining dy 
namo saturation. Regardless of its origin, the linear scaling, 
in any case, nicely explains the shear-rate dependence of the 
dynamo cycle period, which we successfully predict using the 
dispersion relation for the standard a£l dynamo. The striking 
agreement suggests that the positive a effect may be respon¬ 
sible, after all, for the characteristic “butterfly” cycles seen in 
stratified MRI simulations. More supporting evidence for this 
hypothesis comes from the observation that we do see cyclic 
patterns in all studied cases for q, whereas the negative a ef¬ 
fect is found to vanish for q <; 1.5. A definite resolution of 
what effect determines the cycle period, and what effect sets 
the propagation direction of the dynamo pattern will have to 
await quantitative mean-field modeling - taking into account 
the full tensorial structure of both the a and q tensors. 

4.2. Dependence on net-vertical flux 

The inclusion of a significant net-vertical magnetic field af¬ 
fects the amplitude of the turbulent stresses - yet not as dra¬ 
matically as_m_unstratifiedsimulations (compare Fig. [8] with 
Fig. 2 in Pessah et ak||2007 1. For our simulation setup, we 
find a threshold value of - 10 5 , expressed in the mid¬ 


plane plasma parameter. For values smaller than this (i.e., for 
stronger fields), we find an approximate scaling with /?~ l/4 . 
Inclusion of the NVF appears to have a varying effect on the 
various tensor components. We observe some of the coeffi¬ 
cients, such as a xx , and q xy , to be significantly enhanced for 
stronger NVF. The key coefficient a yy entering the oil cycle, 
however, shows the opposite trend, and is weakly quenched 
for strong fields. A similar dichotomy applies to the turbulent 
diffusivity, represented by q yyi which shows a weak quench¬ 
ing near the disk midplane, and is slightly anti-quenched in 
the upper disk layers. We interpret the observed character¬ 
istic to be caused by the faster removal timescale due to the 
stronger wind for higher NVF, such that correlations (caused 
primarily by the Coriolis effect) have shorter time to build¬ 
up in models with stronger NVF and hence a stronger wind. 
The advective helicity flux (Brandenburg et al. |2009 1 will of 
course also scale with the outflow velocity, offering in alter¬ 
native potential explanation for the observed behavior. 

4.3. Scale separation and non-locality 


Motivated by the paper of Brandenburg & Sokoloff ( 2002), 
and with the technological advancement introduced by Bran¬ 


denburg, R(idler, & Schrinner (2008b), we have obtained 


the scale-dependence of the a effect, and turbulent q, for 
stratified, isothermal MRI turbulence with unprecedented fi¬ 
delity. [Brandenburg & Sokoloffj ((2002} for instance found 
q xx - 10 x q yy (see their Figure 12). In contrast, we find 
both components of comparable magnitude. This is likely ex¬ 
plained by their assumption of negligible off-diagonal contri¬ 
butions during the inversion process, which is however not 
warranted in view of our TF results that_show q„ - q xx - q xv . 
Witho ut this restrictive assumption, Brandenburg & Sokoloff 
((2002} find unphysical negative diffusivities for the coefficient 
affecting the radial magnetic field component (see their fig¬ 
ure 10), indicating an inherent problem with their original ap¬ 
proach. 

Using the spectral TF method, for the a yy tensor element, 
we find a wavenumber dependence ~ (\+k 2 /k 2 ) 1 with a char¬ 
acteristic wavenumber of k x r> = 3.7 (in units of 2:t/H) at a 
disk height of z — 2.75 H. The Lorentzian s hape in k space 
has first been ide ntified by Brandenburg et al. (2008b) and Mi- 


tra et al.| ((2009}, and also appears in passive scalar diffusion 
(see Sect. III.C ofjMadarassy & Brandenburg|2010) . 

For the eddy diffusivity q(k), we find a similar dependence, 
but with k^ ’ = 3.3 (in the magnetically dominated corona), 
and 4.8 (near the disk midplane). This illustrates that the 
non-locality in the relation between the mean magnetic field 
and the EMF manifests itself differently in different regions 
of the disk, and this is the first time that both the vertical 
functional dependence and the scale-dependence of these ef¬ 
fects has been derived simultaneously, adding credibility to 
their mutual c onsistency (see the concerns in Brandenbu rg &| 
Sokoloff 2002). With the dominant wavenumber of the mean- 

field dynamo estimated at k, - 0.33, and k c - 3.5, we in¬ 

fer a scale separation close to one order of magnitude, which 

comfortably merits a mean-field description - yet not with¬ 

out concerns that non-local effects may potentially affect the 
properties of the emerging dynamo solution. 

5. CONCLUSIONS 

Its ability to quantitatively explain the period of the ob¬ 
served cyclic field patterns (as a function of shear rate) en¬ 
courages to further explore the aQ mechanism as a key player 
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in the long-term evolution and dynamics of magnetized ac¬ 
cretion disks. The unambiguous determination of the inher¬ 
ent scale-dependence in the dynamo tensors furthermore insti¬ 
gates a new investigation into more realistic non-local mean- 
field models. [Brandenburg & Sokoloff| ( 2002| ) have already 
established that such an approach does affect properties of the 
butterfly diagram in simpler settings. 

A number of issues need to be understood along the way, 
such as the largely unknown (anti-) quenching behavior of the 
mean-field coefficients in the case of a dynamically important 
Lorentz force (even including the possi bility of spontaneous 
symmetry breaking; Chatterjee et al .|201 la' i. A stronger a ef¬ 
fect for stronger mean fields was, for instance, argued to po¬ 
tentially explain the steep scaling of the dynamo cycle period, 
u) cyc , with the rotation fr equency observed in stellar dynamos 
( Chatterjee et al.[201 lb I. 

Alternatively, one can strive to include the effects of the 
Lorentz force explicitly. The latter route would require to ex¬ 
tend the mean-field approach to the momentum equation - 
something that can be regarded as mandatory if one wants to 
include the modeling of turbulent stresses (Qgilvie|2003| 


sah et al.1|2006a|l. Additionally, dynamic al constraints 


Pes- 


rom 


the conservatio n of magnetic helicity (see Hubba rd & Bran^j 
denburg 2012 for a recent account) may govern at least part 


of the dynamics (as witnessed by the dynamically quenched 
mean-field model of |Gressel|2010| l. 

Before global simulations were feasible, a large body of 
work employed the local shearing-box a pproximation (jR 


ley et al[] 1995[ [Brandenburg et al. |1995| to mention two in¬ 
fluential ones) to obtain bulk properties of MRI turbulence. 
In this paper, we also heavily rely on the shearing-box as a 
workhorse for studying differentially rotating turbulence. The 
renewed realization (Suzuki & Inutsuka 2009) of the impor¬ 
tance of outflows a la |Blandford & Paynel (| 1982) challenges 
the usefulness of the local approximation, where some prop¬ 
erties of the outflow are affected by spurious box-size effects 
(Fromang et al,|2013| >. While this may also partly affect our 
results with strong net-vertical fields, we nevertheless argue 
that only combining local and global models can ultimately 
solve the puzzle of understanding turbulent accretion disks. 

All in all, the approach presented in this paper appears to be 
a sensible first step in a program to build a mean-field frame¬ 
work for turbulent flows that will enable to follow the secular 
evolution of global disk dynamics. The fact that the proper¬ 
ties of the turbulent dynamics, as well as the butterfly pattern, 
seen i n local simulations are also ob served in global simula¬ 
tions (Flock et a l.|2012| Jiang et al. 2014| > lends support to this 
approach. 
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